********************************************************************************
********************************************************************************
** Examine how private school students differ from public school students.
** Estimate a multinomial model with three categories:
**	1) Public secondary school enrolment.
**	2) Private secondary school enrolment.
**	3) Dropout.
********************************************************************************
********************************************************************************

preserve

*********
** Examine how private school students differ from public school students
*********
* Load education module from Y3
use "$dataraw_path\nps\2012-2013\HH_SEC_C.dta", clear

keep y3_hhid indidy3 hh_c05 hh_c09 hh_c12 hh_c15 hh_c17 hh_c18 hh_c23 ///
	hh_c28_1 hh_c28_5

gen enrolled = (hh_c05==1) if hh_c05!=.
gen grade = hh_c09
gen private = (!inlist(hh_c12,1,2)) if hh_c12!=.
gen tutoring = (hh_c28_5>0) if hh_c28_5!=.

drop hh_c05 hh_c09 hh_c12 hh_c28_5


* Drop children not in school and without information on school ownership
keep if enrolled==1
drop if private==.


* Keep only students attending lower/upper secondary school
keep if inlist(grade,21,22,23,24)


* Child work, household work, and household education 
merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_E", ///
	keepusing(hh_e08c hh_e32 hh_e50 hh_e64 hh_e69)
drop if _merge==2
drop _merge

egen work_child1 = rowtotal(hh_e32 hh_e50 hh_e64 hh_e69) // student work hours (intensive margin)
gen work_child2 = (work_child1>=1 | hh_e08c==1) // student work dummy (extensive margin)

drop hh_e08c-hh_e69

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_B", ///
	keepusing(hh_b04) nogen // get age of all household members
rename hh_b04 age

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_C", ///
	keepusing(hh_c03 hh_c05 hh_c07) nogen // get school enrolment for all household members

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_E", ///
	keepusing(hh_e08c hh_e32 hh_e50 hh_e64 hh_e69) nogen // get labor information for all household members

bys y3_hhid: egen hhsize = count(age) if age!=.
gen dependent = (age<15 | age>64)
bys y3_hhid: egen dependent_all = sum(dependent)
gen dependency_ratio = dependent_all/hhsize

gen edu_max1 = (hh_c03==2 | inlist(hh_c07,1,2)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // no education
gen edu_max2 = (inlist(hh_c07,11,12,13,14,15,16)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // some primary school
gen edu_max3 = (inlist(hh_c07,17,18,19,20)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // finished primary school
gen edu_max4 = (inlist(hh_c07,21,22,23)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // some secondary school
gen edu_max5 = (hh_c07>=24 & hh_c07!=.) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // finished secondary school or above

bys y3_hhid: egen edu_all_max1 = mean(edu_max1) // share without any education
bys y3_hhid: egen edu_all_max2 = mean(edu_max2) // share with some primary school
bys y3_hhid: egen edu_all_max3 = mean(edu_max3) // share with finished primary school
bys y3_hhid: egen edu_all_max4 = mean(edu_max4) // share with some secondary school
bys y3_hhid: egen edu_all_max5 = mean(edu_max5) // share with finished secondary school or above

egen work_temp1 = rowtotal(hh_e32 hh_e50 hh_e64 hh_e69) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work hours 
replace work_temp1 = 168 if work_temp1!=. & work_temp1>168
gen work_temp2 = (work_temp1>=1 & work_temp1!=. | hh_e08c==1) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work dummy

bys y3_hhid: egen work_all1 = mean(work_temp1) // average work hours in household
bys y3_hhid: egen work_all2 = mean(work_temp2) // share working in household

drop age-hh_e69 dependent dependent_all edu_max1-edu_max5 work_temp1 work_temp2

drop if private==.


* Household consumption
merge m:1 y3_hhid using "$dataraw_path\nps\2012-2013\ConsumptionNPS3", ///
	keepusing(expmR adulteq)
drop if inlist(_merge,1,2)
drop _merge

gen cons_per_cap = expmR/adulteq/1000

drop expmR adulteq


* Present private secondary school correlates
eststo public: estpost summarize tutoring work_child2 hhsize dependency_ratio ///
	edu_all_max* cons_per_cap if private==0, meanonly

eststo private: estpost summarize  tutoring work_child2 hhsize dependency_ratio ///
	edu_all_max* cons_per_cap if private==1, meanonly


* Output of descriptive statistics
esttab public private using "$out_path\tablea2.tex", replace ///
cells("mean(pattern(1 1 1 1 0) fmt(3)) b(pattern(0 0 0 0 1) fmt(3)) se(pattern(0 0 0 0 1) fmt(3))") ///
varlabels(tutoring "Private tutoring" work_child2 "Child working" ///
hhsize "Household size" dependency_ratio "Dependency ratio" ///
edu_all_max1 "Adults without any education" ///
edu_all_max2 "Adults with some primary education" ///
edu_all_max3 "Adults finished primary education" ///
edu_all_max4 "Adults with some secondary education" ///
edu_all_max5 "Adults finished secondary education" ///
cons_per_cap "Consumption per household member")


*********
** Multinomial model
*********
* Load education module from Y3
use "$dataraw_path\nps\2012-2013\HH_SEC_C.dta", clear

keep y3_hhid indidy3 hh_c05 hh_c09 hh_c12 hh_c15 hh_c17 hh_c18 hh_c23 ///
	hh_c28_1 hh_c28_5

gen enrolled = (hh_c05==1) if hh_c05!=.
gen grade = hh_c09
gen private = (!inlist(hh_c12,1,2)) if hh_c12!=.
gen tutoring = (hh_c28_5>0) if hh_c28_5!=.

drop hh_c05 hh_c09 hh_c12 hh_c28_5


* Child work, household work, and household education 
merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_E", ///
	keepusing(hh_e08c hh_e32 hh_e50 hh_e64 hh_e69)
drop if _merge==2
drop _merge

egen work_child1 = rowtotal(hh_e32 hh_e50 hh_e64 hh_e69) // student work hours (intensive margin)
gen work_child2 = (work_child1>=1 | hh_e08c==1) // student work dummy (extensive margin)

drop hh_e08c-hh_e69

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_B", ///
	keepusing(hh_b04) nogen // get age of all household members
rename hh_b04 age

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_C", ///
	keepusing(hh_c03 hh_c05 hh_c07) nogen // get school enrolment for all household members

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_E", ///
	keepusing(hh_e08c hh_e32 hh_e50 hh_e64 hh_e69) nogen // get labor information for all household members

bys y3_hhid: egen hhsize = count(age) if age!=.
gen dependent = (age<15 | age>64)
bys y3_hhid: egen dependent_all = sum(dependent)
gen dependency_ratio = dependent_all/hhsize

gen edu_max1 = (hh_c03==2 | inlist(hh_c07,1,2)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // no education
gen edu_max2 = (inlist(hh_c07,11,12,13,14,15,16)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // some primary school
gen edu_max3 = (inlist(hh_c07,17,18,19,20)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // finished primary school
gen edu_max4 = (inlist(hh_c07,21,22,23)) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // some secondary school
gen edu_max5 = (hh_c07>=24 & hh_c07!=.) if hh_c05!=1 & age>=15 & (hh_c03==2 | hh_c07!=.) // finished secondary school or above

bys y3_hhid: egen edu_all_max1 = mean(edu_max1) // share without any education
bys y3_hhid: egen edu_all_max2 = mean(edu_max2) // share with some primary school
bys y3_hhid: egen edu_all_max3 = mean(edu_max3) // share with finished primary school
bys y3_hhid: egen edu_all_max4 = mean(edu_max4) // share with some secondary school
bys y3_hhid: egen edu_all_max5 = mean(edu_max5) // share with finished secondary school or above

egen work_temp1 = rowtotal(hh_e32 hh_e50 hh_e64 hh_e69) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work hours 
replace work_temp1 = 168 if work_temp1!=. & work_temp1>168
gen work_temp2 = (work_temp1>=1 & work_temp1!=. | hh_e08c==1) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work dummy

bys y3_hhid: egen work_all1 = mean(work_temp1) // average work hours in household
bys y3_hhid: egen work_all2 = mean(work_temp2) // share working in household

drop hh_c03 hh_e08c-hh_e69 dependent dependent_all edu_max1-edu_max5 ///
	work_temp1 work_temp2


* Household consumption
merge m:1 y3_hhid using "$dataraw_path\nps\2012-2013\ConsumptionNPS3", ///
	keepusing(expmR adulteq)
drop if inlist(_merge,1,2)
drop _merge

gen cons_per_cap = expmR/adulteq/1000

gen l_cons_per_cap = log(cons_per_cap)

drop expmR adulteq


* Classify if children are enrolled in public, private, or not enrolled
gen status = 0 if age>=13 & age<=20 & hh_c05==2 & inlist(hh_c07,17,18)
replace status = 1 if age>=13 & age<=20 & inlist(grade,21,22,23,24) & private==0
replace status = 2 if age>=13 & age<=20 & inlist(grade,21,22,23,24) & private==1

label define status_label 0 "Dropout" 1 "Public" 2 "Private" 
label values status status_label  


* Estimate the multinomial logit model
eststo mlogit: mlogit status hhsize dependency_ratio edu_all_max2-edu_all_max5 ///
	l_cons_per_cap, base(2)


* Output
esttab mlogit using "$out_path\tablea3.tex", replace ///
se stats(N r2, fmt(%12.3gc) labels("\(N\)" "\(R^2\)")) compress nomtitles ///
nonumbers unstack noomitted starlevels("" 0.01) substitute(\_ _) b(3) ///
keep(hhsize dependency_ratio edu_all_max* l_cons_per_cap) ///
varl(hhsize "Household size" dependency_ratio "Dependency ratio" ///
edu_all_max2 "Adults with some primary education" ///
edu_all_max3 "Adults finished primary education" ///
edu_all_max4 "Adults with some secondary education" ///
edu_all_max5 "Adults finished secondary education" ///
l_cons_per_cap "Consumption per household member")


restore
